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Abstract 

An anti-noise problem on a finite time interval is solved by minimization of a qua- 
dratic functional on the Hilbert space of square integrable controls. To this end, the one- 
dimensional wave equation with point sources and pointwise reflecting boundary conditions 
is decomposed into a system for the two propagating components of waves. Wellposedness of 
this system is proved for a class of data that includes piecewise linear initial conditions and 
piecewise constant forcing functions. It is shown that for such data the optimal piecewise 
constant control is the solution of a sparse linear system. Methods for its computational 
treatment are presented as well as examples of their applicability. The convergence of dis- 
crete approximations to the general optimization problem is demonstrated by finite element 
methods. 
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1. Introduction 


This article’s optimization problem is motivated by a variant of active noise control in 
acoustics. In physical terms, we consider a wave guide of finite length with the following 
properties: 

(i) the evolution of sound in the guide is appropriately described by plane waves prop- 
agating along its axis without dissipation 

(ii) the sound in the guide originates from initial conditions and sources with small 
support located within the wave guide or at its front ends 

(iii) the system’s total behaviour results from superposition of mutually independent 
sound waves that are generated by the sources and the initial conditions 

(iv) at the front ends of the guide incident waves are partially absorbed and partially 
reflected without alteration of shape. 

The initial conditions at time t = 0 are given and the behaviour of some of the sources 
(offending but uncontrollable sources) is known. Our control problem is to determine the be- 
haviour of the other sources (controllable sources) such as to reduce the noise by destructive 
interference in prespecified regions of the wave guide. 

Most of the above assumptions are used in the literature on anti-noise problems in ducts 
(e.g. [Sw], [TNI, [FW], see also [MI] p.467). Furthermore, in acoustics it is common practice 
to consider harmonic sound fields with time dependences of the form exp(— itot)- In contrast 
to this, the present paper works with non-harmonic sources and fields; for the time domain 
model considered here, a convenient choice for the states and sources are elements in Hilbert 
spaces of real-valued functions. 

The basic time domain model in linear acoustics is the wave equation for the velocity 
potential <f> 

(1.1.1) - c 2 <f> xx {t,x) = /(f,x) + u(t,x), 0 <t <T, 0 < x < L 

(1.1.2) <j>(0,x) = <j>o(x), 4>t(0,x) = tpo(x), 0 <x<L 

with initial conditions <f> o,t/>o • (0, Li) — > R, the constant c > 0 being the speed of sound. As 
is customary in anti-vibration models (e.g. [Sw],[NCEBl],[NCEB2],[P]) we consider point 
sources that are represented by 6 distributions 

n m 

(1.1.3) /(<,*) = u(t,x) = J^Ui(*)^i( x )> 0 < t < T, 

i=2 i=2 

at distinct points 6 (0 ,L) with time dependencies /;,Uj € I/ 2 (0,T). 

According to (iv), the specific acoustic impedances of the boundary surfaces are real 
constants Co>Ci 0- The interaction of a surface and the sound is such that the quotient 
p/vin of the sound pressure p and the incident velocity v in at the surface is equal to pc(i 
([MI] p.259f, p denotes the equilibrium density within the duct). In terms of the velocity 
potential, since p(t, x) = p<f> t (t, x), v(t , x) = - <f> x (t , x), this leads to the boundary conditions 

0<(^i 0) — c^o^i(t) 0) = 9) 

(19! t > 0. 

{ } Mt,L) + chMt,L) = o, 
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We include sources at the boundaries, assuming that sources and reflection combine by 
superposition. For example, the model’s boundary at x — L could be an artificial, non- 
reflecting one (Ci = 1) with independent waves entering from x > L. Let us consider the 
case of an offending source at x = L and a controllable source at x — 0, 

o) - cCo<f> x {t,o) = — ^ - ° «i(t), 

( 1 - 1 - 4 ) * t> 0 

M*,L) + = -^7i(0, 

£ £ 2 (0, T). The factors (1 + £,*)/2c have the effect that the waves generated at the 
boundaries and in the interior have equal amplitudes relative to /*(t),Ui(<). 

p(c 2 <f>x(t,x) 2 + x) 2 )/2c 2 being the acoustic energy density at time t ([MI], (6.2.15)), 

we define the performance index 


(1.3) 


J(ui,...,u m ) 



0 0 


with design parameters r,(f) > 0, i = 1 ,...,m. The weighting function q(t, x) > 0 is 
non- zero on those subintervals of (0, L ) where the acoustic energy is to be reduced. 

The optimal control problem in terms of (1.1) is: Given initial conditions 4>o,ipo and 
functions fi £ L 2 ( 0, T), i = 1, . . . , n, find functions u, € L 2 (0, T), i = 1, . . . , m that minimize 
the functional J wherein <j> is the solution to (1.1). Inserting the solution operators of the 
system into the cost functional J gives a quadratic functional in (uj, . . . ,u m ) whose unique 
minimizer abstractly is represented by use of the adjoint of the operator that maps the control 
functions to the waves they generate. The purpose of the present work is the construction 
of approximate optimization problems that are amenable to numerical computations. 

To this end, equation (1.1) is substituted by a first order system according to the fac- 
torisation dtt — c 2 d xx — (dt — cd z )(dt + cd x ). The components of that system are the two 
components of the velocity potential <f> that are travelling in opposite directions. The delta 
inhomogeneities, generally speaking, involve non-classical but wellknown concepts of solu- 
tion (method of transposition, e.g. [L]). Yet, for the boundary conditions used in this article, 
no directly quotable reference that would cover the wellposedness of the models in section 2 
was found. Instead of following the general theory of [LM], for the one- dimensional systems 
at hand, the existence, regularity and uniqueness of weak solutions (in a sense analoguous 
to [L]) can be proved directly by application of classical methods (d’Alembert’s formula, 
Duhamel’s principle, energy estimate) in the sense of distributions. 

The regularity of the initial conditions and sources in section 2 is taylored precisely to the 
discretized minimization problem of section 3. By restriction to sources that axe piecewise 
constant (in time) and to initial conditions that axe first order splines (in space), the problem 
can be formulated in terms of matrices. This is due to the fact that - according to the simple 
process of propagation and reflection - the value of the cost functional is determined by the 
two components of <f> at the points x } = jcAt at discrete times t* = kAt. Thus, for this type 
of sources and initial conditions, the problem with finite time horizon ( T < oo) is reduced 
to the solution of a finite (possibly large) system of linear equations. 
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The coefficient matrix of this system is sparse. We present numerical methods that avoid 
the zero coefficients. To demonstrate their computational feasability and to illustrate some 
characteristic details of the solutions, four examples are given. 

Finally, twice differentiable initial conditions are interpolated by linear splines and square 
integrable sources are projected onto subspaces of piecewise constant functions. It is shown 
that the corresponding optimal controls converge to the solution of the original problem. 
The proof rests on finite element methods for elliptic problems. 

As is generally the case for wave equations (cf.[DLP],[DW]), the performance of the 
controls is sensitive to time delays. This feature affects the open loop solution of the present 
paper as well as the feedback solution as in [BKSW]. It should be noted that both approaches 
lead to acausal controls in the sense that the behaviour of the offending sources at times later 
than t is used to determine the optimal control action at time t. This is necessary because of 
the global nature of the performance index and the unrestrictive choice of admissible control 
functions. In the examples we can observe the correlation of acausal actions with the finite 
speed of propagation. 

For typical boundary materials that are of practical interest, assumption (iv) is not valid 
([BG],[BPS]). However, modelling and approximating the smearing of incident pulse waves 
by the boundaries (see [MI] p.264f) would destroy the sparsity in the matrix of the discretized 
system and thereby drastically increase the computational requirements. 

On the other hand, we believe that explicitly computable solutions to minimization prob- 
lems for pointwise reflecting boundaries as in (iv) provide detailed insights into the mecha- 
nisms of the optimal control of waves in enclosures. 

Notation. L 2 (a,b;X ) is the Lebesgue space of square integrable functions on the in- 
terval (a, b) with values in X. For real- valued functions, X = R, we simply write L 2 {a,b). 
C([a, b]]X) denotes the space of A-valued functions that are continuous on the closed in- 
terval [a, b] . H l (a, b) is the Sobolev space of absolutely continuous functions h : (a,b) — * IR 
with distributional derivative h' € L 2 (a,b), normed by ||/i||i = (||h|| 2 + ||h , || 2 ) 1 ^ 2 (|| * || being 
the L 2 -norm). As usual, X k = {(aq , . . . , a;*)|aq € X, i = 1, . . . , k}. In a Hilbert space X the 
inner product will be denoted by (•, -)x- We use the abbreviations h(£"^) and h(£ ) for the 
right and left hand limit of the function h at the point £ € R. For an interval / C R we 
write f° r ^e characteristic function of I, i.e. — 1 for t € I, X-f(t) — 0 for t € R \ I- 


2. The first order system 

Any unidirectional linear wave is the sum of two components propagating in positive and 
negative direction. The insertion of the sum of a function of (f — x/c) and a function of 
(t + xjc ) into (1.2) shows, that the boundary conditions (1.2) couple the two components 
by pointwise reflection with reflection coefficients 


-1 < Ri = 


C.--1 

(» + 1 


< 1 , 


* = 0 , 1 . 
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Therefore we consider the system 


( 2 . 1 ) 

( 2 . 2 ) 

(2.3) 


+ c<t>+(t,x) = 0 , 

. 0 < X < L, 0 < t <T; 

4>t (t,x) - c<f> x (t,x ) = 0 , 

0) = Ro<i>-(t,0), L) = R x <j>+{t,L), 0 <t <T; 

<i> + (0,x) = <f>o(x), <t>~(Q,x) = <f>g(x), 0 < x < L. 


This first order system is a generalization of the homogeneous second order equation 
(2.4) <t>u(t,x) - c 2 (f> xx (t,x) = 0, 0 < x < L, 0<t<T 


with initial and boundary conditions (1.1.2), (1.2) in the following sense: If (f) + ,(f>~ are 
functions satisfying (2.1) - (2.3) with initial conditions 


^o"(*) = \m*) ~ Y c J ^ ° d ^ (*) = + y c J ^ ° d & 

0 0 

then <f> = + 4>~ is a solution of (2.4), (1.1.2), (1.2) in the sense of distributions [S2] ; this 

can be seen by the application of the commuting product ( dt + cd x )(d t — cd x ) to (j> + + <f>~ , 
and by verification of the initial and boundary conditions. Classical solutions that are twice 
continuously differentiable evolve only if the initial conditions satisfy certain regularity and 
compatibility requirements. 

2.1. Wellposedness. In context of the present anti-noise problem it is adequate to work 
with continuous initial data, because typical initial conditions axe either silence, or sound 
that has been generated by sources that are driven by L 2 (0, T) functions. Accordingly, we 
define the space of initial conditions 

V = {(4> + ,«r) € H\0,L) 2 \<f>+(0 + ) = R o <t>-(0 + U~(L-) = Rrf+iL-)}, 

endowed with the /T 1 (0, L) 2 -norm. 

Definition 2.1. Given (<f>^, <f>g) £ V, a pair <j>~ ) of functions <f> + , <j>~ : [0, T] x (0, L ) — > 
R is called a solution of system (2.1) - (2.3) if 

(i) <l> + (0,x) = <f>o(x), (f>~( 0,x) = <j>o(x), 0 < x < L; 

(ii) the mapping [0, T] — ► V, defined by t (<j> + (t, •), •)), is in C7([0, Tj; V ); 

(iii) for all x £ (0,L), the functions [0, T] — ► R defined by t ^ <f> + (t,x), t <j>~(t,x) are 
absolutely continuous on [0, T]; 

(iv) for all t £ (0, T), (2.1 ) holds for almost all x £ (0, L). 

Theorem 2.1. Given (<f>Q , <f>Q ) £ V, system (2.1) - (2.3) has a unique solution. With 
£(t) = max{f € N|f < ct/L } the solution is given by 

for £(t) = t odd: 

<l> + (t,x) = Ro^Ro)^ <t>o(-(£ - 1 )L-x + ct ) + (RoR!)^ <}>+((£ + 1 )L + x - ct), 
<j>~(t,x) = Ri(RqRi)~ <$ ((£ + 1)1 - x - ct) + (.Riflo)^ <i>o{-{£ + 1 )L + x + ct); 
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for £(t) = £ even: 

x ) = (R 0 Ri)Uo( £l + *-<*) + Ro(RiRo)$<f>o (~ eL ~ x + c< )> 

<j>~(t,x) = (R 1 R 0 )Ho(-£L + x + ct) + R 1 (RoRi)hU( i + 2 ) L - x - ct )- 

These formulae are understood according to the convention ^o’(r) = 4>q ( r ) = 0 if r € 
R \ (0,L), so that, for every (t,x) G [0, T] x (0,L), at most one of the two terms in each 
formula is non-zero. At the points x = — £L + ct and x = (£ + 1 )L — ct, where the values of 
<j) + and <j>~ are left open by this convention, these values are determined as the left or right 
hand limits, which coincide due to the boundary conditions in V . 

Proof. As to uniqueness, the classical form of an energy estimate (eg.[S2] p.299) is not 
applicable here, because the initial conditions and solutions in Def. 2.1 are not continuously 
differentiable. However, for any solution of (2.1) - (2.3) consider the function 

L 

Eft) = l(<p + {t,x) 2 +<t>-(t,x) 2 )dx. 

0 

Because of (ii), (iii ) , E can be viewed as a distribution whose derivative is given by 

L 

E\t) = j(<)> + {t,x) 2 +cj>-(t,x) 2 ) t dx 
0 

(see [Z], 2.8). Transforming the integrand using (iv) and integrating (<f> + (t, x) 2 + x) 2 ) x 

from 0 to I we get, <p + {t , -) 2 ,<A _ (<, -) 2 bein S absolutely continuous (eg. [HS], (18.16)), 

E\t) = c[(f) + (t, 0 + ) 2 - <f> + {t,L~) 2 + L~) 2 - 0 + ) 2 ] 

= c[{R 2 - l)<j>~(t,0 + ) 2 + (R] - !)<!>- (t,L~) 2 ] < 0. 

This shows that the derivative of E is a non-positive function; therefore ([Sl],Chap.IV) Eft) 
is a decreasing function. Thus, given (<f>Q , (j) 0 ) G V, the difference of two solutions of (2.1)- 
(2.2) is a solution with zero initial conditions i.e. f£(0) = 0. Consequently, for the difference 
of two solutions, Eft) = 0, t € [0,T]. This implies the uniqueness of the solution. 

The formulae for (p + {t, x), x) describe the movement of the two components due to 

propagation and reflection for the (£ + l)th cycle of complete reversion during £Ljc < t < 
(£+1 )Lfc (cf. section 2.2). In light of the convention stated above, it is quite straightforward 
to verify that the given pair (<f> + , <t>~) satisfies (i) - (iv) of Definition 2.1. □ 

With 

H := L 2 (0,T; L 2 (0,L) 2 ) 

we have the following immediate consequence of the formulae in Theorem 1.1. 
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Corollary 2.1. The operator S : V —> H that maps any initial condition in V to the spatial 
derivatives of the corresponding solution of (2.1) - (2.3), S(<f>+ = (<£+(f, ■), •)), 

where <f> + , <f> are given in Theorem 2.1, is a bounded linear operator. 

With regard to sources at the boundaries we consider the homogeneous first order system 
(2.1) with zero initial conditions 

(2-5) <j> + (Q,x) = (f>-(0,x) = 0, 0 <x<L 

and inhomogeneous boundary conditions 

( 2 . 6 ) <P + (t,0+) = R 0 4,-(t,0+), <f>-(t,L-) = R 1( f> + (t,L-) + F(t ) 

where F is related to the boundary source / for the wave equation (2.4) 

^i(t,0)-cCoM^O) = 0, 

(2.7) 1 I { t > 0 


by 


( 2 . 8 ) 



fdr. 


For a source at the left boundary (2.6), (2.7) are to be replaced accordingly. The first order 
system with a point source at £ 6 (0, L) is 


(2.9) 


<t>7 (<> *) - c<t> J(<, x) = cF(t)6i , 


0 < a: < L, 0 < < < T 


where, as in (2.8), 2 cF is the primitive of the point source / in the wave equation 

(2-10) <t>u{t,x) - c 2 (j) xx {t,x) = f(t)6(, 0 < x < L, 0 < t <T. 

The application of ( d t - cd x )(d t + cd x ) to <j>+ + <f>~ , where satisfy (2.9), (2.2), (2.5) 

or (2.1), (2.6), (2.5) shows that these systems are generalizations of (2.10), (1.2) or (2.4), 
(2.7) respectively, with zero initial conditions <f> 0 = V»o = 0 in (1.1.2). 

For £ € (0, L) we denote by the space 


n = {w + ,<n e h'((o,l) \ {(}f i ^ + (o+) = R4f(*> + ),r(L-) = R,r{L-)} 

endowed with the H l norm on (0, L ) \ {£}. 
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Definition 2.2. Given £ € (0,L) and f € L 2 (0,T), a pair (<^ + ,</> - ) of functions <j)+,4>~ : 

[0, T] x (0, L) — ► R with <^> + (0,x) = <f>~(9,x) = 0, 0 < x < L, is called a solution 

a) of system (2.1), (2.6), (2.5) if 

(i) the mapping [0,T] — > defined by t t-» (<j> + (t, •),</»“(<, •)), is in C([0, T]; 

H\0,L) 2 ), 

(ii) for all x € (0, L), the functions [0, T] — * R defined by t <f> + (t, x), t <f) ( t , x) are 
absolutely continuous on [0,T], 

(iii) for all t € (0, T), (2.6) holds, 

(iv) for all t € (0 , T), (2.1) holds for almost all x € (0, L); 

b) of system (2.9), (2.2), (2.5) if 

(i) the mapping [0, T] —*■ V$, defined by t *— ► (<f> + (t, •), <j)~(t, •)), is in C([0, T]; V^), 

(ii) for all x 6 (0,L)\{£), the functions [0,T] — ♦ R defined byt (f>+(t,x), t t— ► <f>~(t,x) 

are absolutely continuous on [0,T], 

(iii) for aii t e (o ,t), <^+(i,e + ) - = r(<,n - r(u + ) = m 

(iv) for all t G (0, T), (2.1) iioids for almost all x € (0, L) \ {£}. 

Theorem 2.2. Let / € L 2 (0, T), £ 6 (0, L). 

a) The system (2.1), (2.6), (2.5) has the unique solution 


oo / 

<f> + (t,x) = Y^R 0 (R,Ro) k F ( 

I — n V 




oo 


— + 1) x[£ — ci,£](— 2fcL — x), 


x) = £(«, Ro) k F ( {2k + ] )L + * + <) x[£ - of, L](-2kL + x). 
fe=o V c / 


b) The unique solution of (2.9), (2.2), (2.5) is given by 


4> + ( t , x) = VffloR, ) k F ( ( — - + f) xK,f + ct](2 kL + x) 

*=o V c / 

+ f; + <) xK - rf,«](-2tl - x), 

fc =0 \ c / 


^-((,x) = y'(R t Rt) k F ( ( 2k . L — + t) xlf - ct.(](-2kL + x) 

fc= 0 \ c / 

+ f; B, (RcR, )*-' F ( { - 2 * L ± f + () xlf, f + ctl(2fct - x). 
fc=i V c / 

Proof. The formulae are constructed by folding the solution of the problem on all of R (no 
boundaries) into the domain (0,L) with reflection coefficients Ro,Ri. Note that for any 
t > 0 all the sums are finite, because for large enough k € N all characteristic functions are 
evaluated outside their support. Thus, the solution-properties of the given pairs ,<j>~) 
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can be verified by considering the sums term by term. F being absolutely continuous, 
the regularity with respect to x and t follows from the fact that F is zero wherever the 
characteristic functions have time dependent steps. That the norm of (<f> + (t, •), •)) 

varies continuously with t is due to the steady propagation with finite speed c (cf. section 2.2). 
As to uniqueness, the difference of two solutions of one of the inhomogeneous problems 
(extended continuously at £ in case b)) is a solution of the homogeneous problem (2.1), 
(2.2), (2.5) and thus equals zero by Theorem 2.1. □ 

For £ € (0, L) let denote the operator that maps any inhomogeneity / € T 2 (0, T) to 
the restriction to (0,L) \ {£} of the spatial derivatives of the solution of (2.9), (2.2), (2.5), 
i.e. 

S(f(t) = •)) € L 2 ((0, L) \ (0) 2 , 

where fa , fa are given in Theorem 2.2 b). The point £ is spared out to exclude <*>- terms in the 
derivatives of fa(t , •), fa(t y •), whereas the velocity potential fa(t, •) + fa(t, •) is continuous 
at £. Throughout we can choose = (0,0) as a possible extension to (0, L). For 

notational convenience we let = I, rji =0 and denote by S 0 and S L the operators that map 
boundary sources / € T 2 (0, T) to the spatial derivatives (t, •), </> J (t, •)) of the solution of 
(2.1), (2.6), (2.5), where (for the right hand boundary) fa, <j> are given in Theorem 2.2 a). 
An immediate consequence of Theorem 2.2 is 

Corollary 2.2* For any £ G [0, L], : L 2 (0 y T) — > H is a bounded linear operator. 

2.2. Propagation and reflection. The components of the solutions given in section 2.1 

evolve by shifts with speed c in positive and negative direction combined with inversion and 
reduction by i?o, R\ the boundaries. The source / in (2.10) generates a wave that is 
symmetric about £ G (0, L). More specifically, let the derivatives of the two components 
of the velocity potential (z/>+(t, •), 0“(t, •)) := be given at some time t > 0 and let 

0 < r < T£ = min (£,£ — f)/c. Then, at time t + r 

fa(t + r,x) = + t + r ) x[£,f + cr](x) 

^ f fa(t,x — cr), ct < x < L 
/ 211 \ \ -Rotp-(t,cT - x), 0 <X<CT 

i>~{t + r > x) = 2^2 f + t + t\ x[f - cr, £](x) 

^ f fa(t,x + cr), 0 < x < L — ct 

\ —R\fa(t, 2 L — x — cr), L — ct < x < L. 

This is to be understood almost everywhere in (0, L ) \ £. For a boundary source at £ = 
0(£ = L ) the first term in fa(t + r, x) (fa(t + r, a;)) is to be deleted with tq = ti = Ljc. 
For the pair (fa(t, •), fa(t, •)) := S(fa , (2.11) applies for 0 < r < Ljc with / = 0. 

All this can be checked by differentiating and regrouping the formulae of section 2.1 and 
strict adherence to the convention at Theorem 2.1. 

2.3. The minimization problem. For superposition of solutions to initial conditions 
<t> o > fa and multiple offending sources /, ,...,/„ at the points 6 , . . . , £„ € (0, L) define g G H 
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by 


?« = sw3\«<) + £>«.*(<)• 

1=1 

For the controllable sources located at r/i , . . . , r) m G [0, T) we write 

u = (u u ...,u m )€U:=L 2 (0,Tr 

and denote by B the bounded linear operator U — ► H 

m 

Bu(t) = Y^, s vM t )- 

i=l 

In this notation the spatial derivative of the solution of the first order system that corre- 
sponds to (1.1) is g 4- Bu. Rewriting the functional (1.3), note that for solutions {4> + , 4>~) 
of (2.1) 

c 2 {<j> + + <f>~) x (t,x) 2 +(<t> + + <t>~)t(t,x) 2 =2 c 2 [<f>+(t,x) 2 + <j>~(t,x) 2 ], 

for all t G (0, T) and almost all x G (0, L). Therefore, we are considering the following 
minimization problem: Given (<f>o , 4 >q ) G V and fi, . ■ ■ , f n £ L 2 (0 , T), determine u G U 
such that 

(2.12) J(u) = min{J(u) : u G U), 
where 

(2.13) J(u ) = (g + Bu, Q(g + Bu))h + («, Ru)u- 

Here Q and R are bounded selfadjoint operators on H and U resp., defined by 

Q(<f> + ,<f)-)(t,x) = p(q + (t,x)<j> + (t,x),q-(t,x)<l>~(t,x)) 

R(u I , . . . , u m )(t) = (rj (t)uj (t), . . . , r m (t)u m (t)). 

The two components of the waves may be weighted separately, but the functionals in (1.3) 
and (2.13) coincide if q + = q~ = q. We make the assumption that the weights q + ,q~ G 
L°°(0,T; L°°(0, L)) and G L°°(0,T) are such that 

(2.14) R + B*QB> 0, 

where B* : H — > U is the adjoint operator of B. Then the unique minimizer of the quadratic 
functional J is given by (eg. [B], (5.2.4)) 

u = -(R + B*QB)~ l B*Qg, 
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3. Discretization 


In order to obtain a discrete version of the minimization problem, we consider a grid of 
characteristics of (2.1) in the (<,x)-plane with knots at the points ( tk,Xj ), where xj = jAx , 
j = 0, . . . , N and <* = kAt, k = 0, . . . , K. We assume L = N Ax, T = K At and Ax — cAt. 
Furthermore all sources should be located at meshpoints, i.e. {£i, • • • , £n,*/i, • • • ^m} C 
{x 0 , . . . ,x/v}. 


3.1. Piecewise constant sources and spline initial conditions. Define the finite- 
dimensional spaces 


Uk = {f ■ [0,T] -> R | / is constant on [<*- 1 , 4 ), k = 1 

Vn = <t>~) G V | (f> + , <f>~ are first order sphnes with respect to {xj}y[L 0 }, 

Hn = {(V ,+ , */>-) € H | xf + ,xp~ are constant on [xj_i , Xj),j = 1, . . . , N}. 

A first order spline with respect to {xj}j)L 0 is a continuous function [0, L] — * R that is a 
polynomial of degree one on [xj -\ , Xj], j = 1, . . . , N. Being subspaces of L 2 ( 0, T) and H, Uk 
and Hn are endowed with the corresponding L 2 -products. For the kth basis element in Uk 
we choose x[<*-i, <*), k = 1 ,K and we use the pairs (x[^j-i , Xj ), 0) and (0, xKr-i , x j)) 
for the jfth and ( j + lV)th basis elements in Hn , j = 1, . . . , N. 

Lemma 3.1. If (<f>Q, <f>$) E Vn and /* E Uk, i = 1, . . . ,n then g(tk) 6 Hn, k = 0, . . . ,K. 
If Ui E Uk, i = 1, • • • ,m, then Bu(tk ) € Hn, k = 0, . . . ,K. 

Proof. Because of Ax = cAt, the statements of the lemma follow from (2.11). □ 

For initial conditions (4>q,4>o) G Vn and sources /, £ Uk the coordinate vector g of 
the ( K + l)-tupel (<7(<o), • • • ,g(tK)) G H has 2N(K + 1) real entries. The coordinate 
vector u of u = («i E Ujf is in R mft \ The operator that maps any u E U™ to 
(£u(< 0 ), . . . , Bu(t K )) has a 2N(K + 1) x mK matrix representation S. With this notation, 
the time-sampled spatial derivative of the solution to the first order system that corresponds 
to (1.1) is represented by the vector j + Bu. 


3.2. Discrete minimization. We replace the minimization problem (2.12) by the 
following one: given (4>q , <f>Q ) E Vn and /i , . . . , /„ E Uk, determine u E U™ such that 


(3.1) 


J(u) = min{ J(u) : u E U%}, 


where J(u ) again is defined by (2.13). 
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With the numbers 


pi 


= » / 


t k Xj-l+c{t — tk-l) 


i 


q + (t, x)dxdt , 


Pk 


j q + (t,x)dxdt, qt,j=P j 

ffc-1 Xj 

_ tk x i 

j = p j J q~(t,x)dxdt, qkj=P J J q~(t,x)dxdt, 


tic -1 x j - 1 

tk Xj~c(t-t k - l) 


tk - 1 x j - i 
tfc 


— c( t — t — l ) 


ri.fc = J ri(t)dt. 


tk - 1 


we define the 2 iV x 2 N matrices 


Qo = diag( 9 fu • • • » 9 *jv>Pi,i» • • • ’Pi,n)' 

Qk =diag + 9 if+i l i*-”»P]tw + 9 lf+i,N»Pib+i,i + 9 fc,i’ 

QiC = diag(p^j , • • • ,P/f,An 9jf,i » • • • > )• 


fc = 1, . . . , A" — 1, 


and the square matrices 


Q = diag(Q 0 ,- •• ,Qk), 

R = diag(ri,i,. . . , . . . ,r m ,i , - • • ,r m ,K)- 


of size 2 N(K + 1) resp. mK. 

A minimization problem written in terms of these matrices is: determine u € lR m such 
that 

(3.2) J(u ) = min{ J(u) : u 6 

where 

J(u) = (fif + Bu) T Q(g + B_u ) + u T Ru. 

We assume f? + B^ QB_ > 0. 

Theorem 3.1. The coordinate vector u of the solution u of (3.1) is the solution of (3.2). 

Proof. Let u € Uf) and u = c.ol(ui ; i , . . . , iti./c, • • • > ^m,i , • • • > Um,K ) be its coordinate vector. 
Then 

m K 

( u,Ru)u = 

i=l k= 1 

Given spline initial conditions and piecewise constant sources, let h(t) = (h + (t, •), h (t, •)) = 
(g + Bu)(t) and denote the coordinate vector of (h(t 0 ), . . . , h(tx)) £ #n + 1 by A = c°l(h^j , 
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• • • > ^o,ni ^o,n • ■ • 5 ■ • • i ^K t 1 5 • • • > ^k,N'> ■ • > h K N ). The movement of the two com- 

ponents of /i(t) over (xj-i, Xj) during t*-i <t <tk is indicated in Fig. 3.1. 


Xj - c(t - t k _i) 


x j-l +c(t - <*_! ) 


h h 


k-l,i 


h t , 






Fig. 3.1: Piecewise constant components at time t G (t/t-i , tjt)- 


We see that 




J j (q+(t,x)h + (t, x) 2 -(- q ( t,x)h (t,x) 2 )dxdt 

*k-l x ) - 1 

= + q t, h t 2 h} + PliK- i,i + iK 2 

for all k = 1, . . . , K and j = 1 , . . . , JV. Using this, we get 


K N 


t* x i 


(h,Qh) H = f j (9 + (<> a: )J l+ (*» a 0 2 + 9 ( t,x)h (t,x) 2 )dxdt 

1 1 i -J J 


k = \ >= 1 


TV 


tk - 1 x j - 1 


TV 


~ E q hj h 0J + Pi,j h 0 2 j + ^2PK,j h K 2 j + QK,j h K 2 j 

i=l j=l 

+ E ErfX; + + EEW-u 

fc=l >=1 Jt=2>=l 

= A t QA. 

Therefore, J(u) = J(u). □ 

Thus, the vector representation u of the solution to (3.1) solves the linear system 


(3.3) 


(R + B t QB)u = -B t Q g. 
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4. Numerics 


Without reference to the particular structure of system (3.3), its solution is a standard 
task. However, since the discretizations in time and space are coupled, A t = Ax / c, the size 
of system (3.3) is large if high spatial resolution is required and/or if T is large. Then the 
number of entries in the full matrix B_ is beyond the memory capacity of standard hardware. 
To overcome this problem, we present efficient algorithms for the computation of the right 
hand side of (3.3) and the nonzero entries of the sparse matrix B QB. 

The vector g can be produced by discrete simulation of the propagation and reflection of 
the components 4 > 'x ■ <7(1 • 27V) is given by the initial conditions. Instead of shifting the 
arrays that hold at each time step, it is more efficient, in particular when only a few 

of the time steps have to be assigned (Example 4.3), to keep the arrays, except of replacing 
the entry for </>+( 0 + ) by -Rq times the previous entry for </>“(0 + ) and the entry for <f> x (L ) 
by -Ri times the previous entry for To accomplish this, an integer variable for the 

current index of <f>t(0 + ) is initialized to 1, decreased by 1 at each time step, or reset to TV if 
it is 1. The same boundary index is applicable for both components, if the array for (p x is 
arranged in reversed order. To keep track of the indices where to add the contributions of 
the offending sources /,, we have to initialize and update integers that point to the source 
locations i = 1, . . . , n. The assignment of the entries of <j >+ , <j>~ to g(2kN + 1 : 2(fc + 1 )N) 
gives the fcth group of g , k — 1 , . . . , K . For the subsequent multiplication with B_ is it 
convenient to store the <f>~ groups in reversed order. Q being a diagonal matrix, Q g is 
obtained from g by 2 N(K + 1) multiplications. 

For the computation of —Bj'Q g, consider the structure of B_ . Using 2c as unit for 
the source amplitudes t<i(f), we omit the factors i n (2.11). Then the A;th basis 

elements of U K produces, at time t k , a </>" pulse to the left and a <j>t pulse 

— x[xj, Xj+i ) to the right of the source location (boundary sources produce only one such 
pulse). These pulses travel one space index per time step and are reflected at the boundaries. 
Thus, the jfcth row of a source block in starts with k 2iV-groups of zeros and continues 
with K - (k - 1) 27V-groups each containing (at most) two nonzero propagating entries. 
Initializing these entries to 4*1 according to the source location, and updating their indices 
and amplitudes according to propagation and reflection, the fcth entry of the ith source in 
-B T Qg is accumulated by adding the products of the pulse amplitudes with those entries 
of -Q~g that are determined by the updated indices. Repeating the entire procedure for 
i = 1, . . . , m, -B T Q g is computed avoiding the zeros in B T (which is a mK x 2N(K + 1) 
matrix with at most mK x 2 K) nonzeros). 

Next, we describe an algorithm to compute Bj' QB_ for time independent weights q (t,x) 
= JKiEi Then Q = diag^Qi, Qi, . . . , Q, , |Qi) with Q 1 = diag ( 9 +, 

e R2N * The al g° rithm is based on the fact that two waves emer S in S 

from two sources contribute to an entry in B only if the pulses that are generated by 
the basis elements of Uk overlap during their travel forth and back in the duct. Two pulses 
overlap starting with the generation of the later one, or never. Therefore B J3 consists of 
K x K blocks that are banded, each nonzero diagonal containing the accumulation of the 
products of two overlapping travelling pulses. 
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An entry in a diagonal can be computed from the one south east to it by adding the 
contribution that occurs because both sources are active one cycle earlier. For doing this, 
it is convenient to imagine the overlapping pulses travelling backward in time, adding con- 
tributions to what has happened later. At each time step, the product of the overlapping 
pulses has to be multiplied by q* or 5”, according to the current direction of propagation 
and the current location of the overlapping pulses (by half the value of q+ or qj for the last 
time step). For the four lower diagonals near (resp. on) the main diagonal of the block, 
the location of the entry in the most southern row is determined by the distance of the two 
sources of the block and their distances to the boundaries; the product of the overlapping 
pulses is 1, R 0 , R \ , RqR\ respectively. These four diagonals are filled up one after the 
other: The elements to the north west of the most southern entry are computed until the 
overlapping pulses arrive at the boundary. From then on the diagonal entries are computed 
in groups of N elements corresponding to the propagation back and forth in the duct. At 
the beginning of each group the product of the pulses is multiplied by the square of the 
appropriate reflection coefficient. The diagonal is complete when the most western colu mn 
of the block is reached. 

Each time a nonzero entry to one of these four diagonals is determined, the entries within 
the block in the same row and 2 j N columns to the west of it are computed by multiplication 
with (R 0 R 1 ) 3 . These represent the interaction of a source that was active at the same 
location but j full periods (= j2N cycles) earlier than the one just considered - its pulse 
was reflected 2 j times before overlapping with the later pulse. 

Analogously, the diagonals above the main diagonal of the block are built up from south 
east to north west with the full periods north of the current entry. 

These methods, coded in Fortran 90 (NAGWare compiler 1.2), create disk files that 
contain the right hand side and the nonzeros of the lower triangle of the coefficient matrix 
in (3.3). The latter one consists of lines of the form row index, column index, value of 
the entry. The files are loaded into Matlab 4.1 where the coefficient matrix is set up as a 
sparse symmetric matrix A and the right hand side is a full vector rhs. The solution u is 
then obtained by the Matlab command u = A \ rhs that invokes sparse matrix arithmetics 
[M]. The generation of the optimal waves is again done in Fortran 90, analogously to the 
generation of g. Finally, we use Matlab for the graphics. 

In the following examples, t is given in seconds, x in meters, c = 344. The relation of 
the diagonal elements of Q, denoted by q ± (x), to the weights q ± (t,x) is given above. We 
use time independent functions for r,(<), so that the entries of R are r,- = r,(t)At. The step 
functions /,•(<), «,•(<) are measured in units of 2c 2 . Then J(u ) = Jju) = h T Qh + u T R u. 

Example 4.1. L = 10, R 0 = l, Ri = 0.5, N = 40, K = 137, n = 1,6 = 4,m = l,r/i = 
0, rj = 0, q±(x) = 1, zero initial conditions. This introductory example is comparable to 
the problems in [CNE]. The offending source f\ = x[ti9, £20) emits a positive and a negative 
pulse at time <20 = 20 Ax / c ~ 0.0145. To reduce the rightgoing pulse of the offending 
wave, the optimal control source acausaly emits a pulse wave at time 6: ui(<) = —0.68 for 
t G [£3^4) • T° annihilate the pulse that comes from in negative direction, iii(t) = —1 
for t G [^35, t 36 ); to annihilate the reduced pulse that is reflected at the right boundary, 
^1(6 = —0.16 for t G [t83 ? f84)‘ After that <j> is constant, silence. The 3D graphics show the 
velocity potential at the grid points ( tk,Xj ) connected by straight lines that are parallel to 
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the axes. In the plane of the t , x axes a contour plot indicates curves of constant <f>(t,x). 
Compared to the uncontrolled process J( 0) = 116.22, the optimal acoustic (r i = 0) energy 
is J(u) = 26.88. 

Example 4.2. The configuration is as in the previous example, except that the control 
source now is located in the interior, 771 = 2. The optimal control thus generates leftgoing 
pulses so that the annihilation of <f>± is never completed. J(u) = 56.95. 

Example 4.3. L = 10, Ro = f?i = 0.2, N = 25, K = 148, n = 1, £i = 10, m = 2,771 = 
3.2, n = 0.5,772 = 7.2, r 2 = 0.5, ^(i) = x[0, 3.2). The offending source at the right bound- 
ary runs at 50 Hz, fi(t) = sinlOOTr** for t € [t*-i,<*)- As harmonic initial conditions we 
take from the last time step of the (uncontrolled) generation of g that starts form 

pre-initial conditions <£+( — 10, •) = <£j( — 10, •) = 0 with f\ running for —10 < t < 0. The 
optimal control causes a transient (t > 0) toward another periodic state with a flat velocity 
potential on x (0,3.2). Note that, because of r 2 / 0, u 2 switches to 0 as soon as its 

signals cannot influence the waves in (0, 3.2) any more. We get .7(0) = 659.46, J(u) = 20.65. 
The flatness of <f> in (0, 3.2) can be influenced by the magnitudes of r, , r 2 relative to q ± . For 
our choice of r 1 , r 2 we see some variation of <j> throughout. For smaller design parameters 
ri,r 2 the oscillations of <j> become invisible on (fis,**') x (0,3.2). However, r 2 = 0 would 
yield a singular coefficient matrix since q ± (x) = 0 on (3.2, 10). 

Example 4.4. To test the applicability of our methods to somewhat more realistic acoustic 
events, we sample the first 2 seconds of a wellknown theme of classical music, fi(t) = 
£fc=i x[|> ^ tl )sin3927rt + x[l,2)sin3127rt, t 6 [0,2], and feed it to a source at £1 = 2.5 
in a duct of length L = 5 with reflection coefficients Ro = 0, R\ = 0.7 and zero initial 
conditions. We seek to reduce the sound near the boundaries: q ± (x) = x(0, 1) + x(4, 5). Let 
two control sources be located at 771 = 1, r/ 2 = 4 weighted by r 1 = r 2 =0. We would like to 
have spatial resolution of about 30 grid points for the shortest wave length of the offending 
wave, 344/196/30 = 0.059, so we choose Ax = 0.05, = 100. Then T = 2 = KAt leads to 

K = 2c/ Ax = 13760. 

The optimal control annihilates the components of the wave that cross 771 , t] 2 in outgoing 
direction. This amounts to total reflection at the boundaries of (771,772), which, roughly 
speaking, in general leads to increasing sound in its interior where the offending source 
is located. This is not the case for the symmetric configuration chosen here. We give 
graphics for the transients at t — | and t = 7. At the end of the first tone, the graphics 
show a transient toward a smoother velocity potential. At the end of the third tone, the 
controls generate waves that are absorbed by the offending source, so that <f> = 0 (although 
n = r 2 = 0 and q ± = 0 in (771 , 772)) until the fourth tone of lower frequency begins at t = 1 
(the contour curves reveal computational inaccuracies). Because of the totally reflecting 
behaviour of the control sources, there is no damping in [771,772] and the process does not 
converge to a periodic state as t — >2. In fact, the choice of Ro, R\ does not affect the optimal 
control for this configuration; however, R 0 = 0 drastically reduces the number of nonzeros 
in (3.3). The full matrix 13 would require over 600 gigabytes of storage (8 bytes per entry), 
whereas the sparse matrix A of density 0.0003 has 2.4 megabytes (including the memory for 
the indices). The sparse arithmetics to solve Au = rhs require less than 10.3 megaflops. 
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5. Approximation and Convergence 

We again, assume that L, T and the location of the sources are such that there exist integers 
N,K £ N with Q {jL/N,j = — 1} and K = cTN/L. 

For £ £ N and <f> £ //^(O, L) let st<j> be the first order spline interpolation of <j> with respect 
to the mesh {jL/£N}*= 0 , i.e. 

s t <f>(jL/£N ) = <j>(jL/£N), j = 0,...,£N. 

There exists a Co > 0 such that for all £ £ N and <j> £ H 2 (Q,L) 

11$*^ — ^ii//i(o,D < 

(e.g. [Sch], Th. 2.5). Moreover, let 717 denote the orthogonal projection of X 2 (0, T) onto 
Uik , i.e. 

kT/tK 

tK f 

= y J fdr for *€[(*- 1)T/£K, kT/£K), k = l,...JK. 

(k—l)T/tK 

For aU / € L 2 (0,T) (cf. [BB], p. 176) 

(5.2) |K// — /||l*(o,t) — ► 0 as £ — ► 00 . 

Given /i,...,/ n € L 2 (0,T) and sufficiently smooth <f>Q,<f>Q, the solution u of (2.12) is 
approximated by the functions ut, that are the solutions of (3.1) wherein (se<$,st<f>Q) and 
tt//i , . . • , 7T// n are used for initial conditions and sources and the minimum is taken over 
Uff(. More precisely we have 

Theorem 5.1. If (<f>Q,<f> 0 ) € V fl H 2 ( 0, 2) 2 or if (<f>Q,<f>o) £ Vi 0 n for some £0 € N then 
||u — u<||t/ — > 0 a s £ —* 00. 

n 

Proof. With gt = S(st<pQ ,se<f> 0 ) + ^2 the optimal control ue £ is characterized 

»=1 

by the variational equation 

a(ut, u ) = bt(u) for all u £ Uff. 

where a : U x U — ► R, 6 e : U — » R are symmetric bilinear, resp. linear, continuous forms 
given by 


o(u, v) = (u, (R + B*QB)u)u 

bt(u ) = -(B*Qgt,u)u. 

This is so, because u t minimizes the functional J t (u ) = a(it, u) - 2 bf(u) + (g t , Qg e ) H in U 71 
(cf. [C], Theorem 1.2; the constant term is included here for consistency with (2.13), but it 
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does not affect the minimization problem). Similarly, u minimizes J(u ) = a(u,u) — 2b(u) + 
(d,Qg)H in U with 

b(u) = -(B*Qg,u)u 


and thus it is characterized by 


a(u, u ) = b(u) for all u £ [/. 

By assumption (2.14) a is [/-elliptic and we can apply the first Strang lemma ([C], Th. 26.1): 
there exist constants ci , c 2 that axe independent of £, such that 

||u - u t \\u < ci inf{||u - u||[/ : u 6 U™ K ) + c 2 sup{|6(u) - 6 < (u)|/||u|| [ / : u € U™ K } 

The first term on the right hand side tends to zero by (5.2). The second term is estimated 
using Schwarz’ inequality and the boundedness of B*Q 

I (B*Q(g - gt ), «)i/|/IMIt/ < ll#*< 2 llll <7 ~ 9t\\H 


where 


II0 - gt\\H < ll-S'HIK^o » ^0 ) - 0 )lk + ^2 - **/sIU*(o,t) 

i=l 

so that (5.1), (5.2) imply the result. In case (^, <f>Q ) € Vt 0 n we have st<j > ^ = <f>^ , st<f> 0 = <f> 0 
for l > Iq. □ 

The second case in Theorem 5.1 is applicable when the problem is broken up into pieces 
over subintervals of [0,T]. But this is not equivalent to minimization over [0, T] at once. 

How do the discrete controls ui perform when applied to the original data (<7^,<^ 0 ) € 
V n H 2 (0,L) 2 and sources fi e T 2 (0, T)? The answer (in the sense of convergence) is 
contained in the theorem: J(ui) —* J(u), since J :U —* R is continuous. 

5. Conclusions 

By decomposition of the one-dimensional wave equation with point sources and pointwise 
reflecting boundary conditions and by appropriate discretization in time and space, the 
minimization problem with piecewise constant sources was reduced to the solution of a 
sparse linear system. It was shown that the waves that are generated by discrete forcing 
functions and initial conditions are sufficiently regular so that the problem is wellposed. Four 
examples demonstrated the applicability of algorithms that efficiently set up the sparse linear 
system; such examples provide insights into the optimal control of waves in ducts. Finally, 
it was proved that the solutions of finite dimensional discretizations converge to the solution 
of the minimization problem with square integrable sources. 
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